Synchronization transition in scale-free networks: Clusters of synchrony 
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We study the synchronization transition in scale-free networks that display power-law asymptotic behaviors 
in their degree distributions. The critical coupling strength and the order-parameter critical exponent derived by 
the mean field approach depend on the degree exponent X, which implies a close connection between structural 
organization and the emergence of dynamical order in complex systems. We also derive the finite-size scaling 
behavior of the order parameter finding that the giant cluster of synchronized nodes is formed in different ways 
between scale-free networks with 2 < X < 3 and those with X > 3. 
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I. INTRODUCTION 

Much attention has recently been paid to integrative ap- 
proaches to complex systems in diverse fields. Graph- 
theoretic analyses on social, biological, and informational 
networks consisting of relevant elements (nodes) and inter- 
actions (edges) have found common features such as short 
average distance between nodes, high clustering, and scale- 
free (SF) topology, which means that the degree distribution 
is of a power-law form lily]. These structural characteris- 
tics have been found to be evidence of competitive evolution 
processes 0] as well as ground for the robustness of com- 
plex systems against external attacks or internal errors jj. 
However, the relation of the structural organization to dy- 
namical complexity remains to be explored. Recent works 
on epidemic spreading 0], avalanche dynamics 0, 01, and 
reaction-diffusion processes 1 8, 9] in complex networks iden- 
tified the significance of network topology in dynamic behav- 
ior. Here we focus on collective synchronization phenomena 
that appear in various physical and biological systems includ- 
ing Josephson junction arrays, cardiac pacemaker cells, flash- 
ing fireflies, and brains 1 10] as well as are required in con- 
trolling artificial systems |11]. While each element of these 
systems can be described simply by a limit-cycle oscillator, 
the coupling topologies are not always regular. For instance, 
the neuronal network with syna ptic coupling in human brains 
turns out to have SF topology il2l fl3ll . whose synchrony is 
essential in the recognition of spatio-temporal patterns (3. 
Also in parallel simulations where tasks are distributed among 
many processing elements, weak random couplings among 
them can prevent the spread of virtual times from diverg- 
ing llltl . Our interest is in how the emergence of dynamical 
order is intertwined with structural complexity. 

The linear stability analysis of the synchronized state in 
complex networks of oscillators has shown that the synchro- 
nizability, quantified by the ratio of the largest and the sec- 
ond smallest eigenvalues of the Laplacian matrix, would be 
suppressed by a few "center" nodes being overloaded by the 
traffic of communication 1151 Ma H7I1 . However, the ef- 
fects of topological features on the phase transition from the 
desynchronized state to the synchronized state demand fur- 
ther investigation. In small-world (SW) networks (JJ] , the 
synchronization transition is observed at a finite critical cou- 
pling strength and associated critical exponents are shown to 



be equal to those of the globally-coupled case fl8l [T^ EbTl 
meaning both belong to the same universality class 12111 . In 
this paper we report on the synchronization transition in SF 
networks of limit-cycle oscillators. The degree distribution 
of the SF networks takes a power-law form p c i(k) ~ k~^ for 
k 3> 1 • We find by the mean-field approach that the critical 
coupling strength is much smaller than those of completely 
random graphs 1 22] or S W networks 1 1 ] that have Poisson de- 
gree distributions, and the order-parameter critical exponent 
varies continuously with the degree exponent X as long as 
2 < X < 5, illustrating the relation between the structural or- 
ganization and the emergence of dynamical order. We also 
show how the mean-field approach enables to compute the or- 
der parameter in the critical regime for finite system size and 
derive the finite-size scaling behavior of the order parameter 
finding different nature of the synchronization transition be- 
tween SF networks with 2 < X < 3 and X > 3. Our numerical 
data confirm the analytic results. 

The paper is organized as follows. The model for the cou- 
pled limit-cycle oscillators on a network is introduced and the 
order parameter is defined in the following section. In Sec. Mil 
we use the mean-field approach to set up self-consistent equa- 
tions and find the critical point and order-parameter criti- 
cal exponent. The finite-size scaling behavior is derived in 
Sec. lIVI We summarize and discuss our work in SecfVl 



II. MODEL AND ORDER PARAMETER 

We consider N coupled limit-cycle oscillators whose phases 
{(|),(f)|/ = 1,2,... ,N} evolve according to 



dt 



<?;;Sin 



(<l>i-f/), 



(1) 



where J is the coupling strength, ey = 1 (0) if oscillators i and 
j are connected (disconnected), and (k) denotes the average 
degree £,-&,-/iV with k[ = Y*j e ij- The natural frequencies co,-'s 
are distributed following the Gaussian probability distribution 
g{(Hi) = (2z)~ l / 2 e~ (l> l/ 2 . A system of globally-coupled oscil- 
lators evolving by Eq. Q, i.e., the case of e, 7 = 1 for all / ^ j 
and (k) = N — 1, defines the Kuramoto model I18il and has 
been used as an exactly solvable model of collective synchro- 
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FIG. 1 : Order parameter r = lim,^oo r{f) versus the coupling strength 
J for N = 1600, (k) = 4, and X = 6.4, 3.6, and 2.4. Inset: Time evo- 
lution of r(t) for X = 3.6 with J = 0.0,0.4,0.8, 1.2, 1.6,2.0,2.4,2.8 
and 3.2 from bottom to top. 



nization ficll . Here we are interested in sparse SF coupling 
networks with (k) = O (1) and X > 2. 

Synchronization may be quantified by the amplitude order 
parameter r(t) defined by 



r(f) 



1 

N 



■<M0 



The average phase 8(f) satisfies dQ(t)/dt = CO with CO — 
(l/N) f° r a given frequency realization {co,}. To com- 

pute r(f), we performed a numerical integration of Eq. Q on 
SF networks that are generated using the static model intro- 
duced in Ref. l23ll . r(t) is averaged over several hundreds of 
different frequency and network realizations, {co,} and {e,y}. 
As shown in the inset of Fig.Q r {f) saturates to a value r in 
the long-time limit. We plot r for N = 1600 as a function of 
J in Fig.[l] evaluated by r = (2/T)Y^ =T j 2+i r (0> where we 
set T so as to guarantee that the saturation of r has already 
begun at time T /2. The numerical results imply that the tran- 
sition from the desynchronized state (r = 0) to the synchro- 
nized state (r non-zero) in the thermodynamic limit N — > °° 
exhibits a dependence on the coupling topology characterized 
by the degree exponent. The critical phenomena related to the 
synchronization transition in SF networks can be understood 
by the mean-field approach. 



III. SELF-CONSISTENT EQUATIONS AND CRITICAL 
EXPONENT 

In this section, we obtain the mean-field solution to Eq. Q 
and find the critical point and order-parameter critical expo- 
nent. Assuming the same magnitude of effective field be- 
tween each pair of coupled oscillators, one obtains the one- 
body differential equation from Eq. Q as d§i/dt = CO, — 



(// (k))kir(t) sin (ft;, where the effective field r(t ) satisfies 
lf = iI*Li^exp(^) 



Li=i Lj=i Cij 



(3) 



and the relation 9(f) — > in the thermodynamic limit is used. 
?{t) will be computed in a self-consistent way. 

In the stationary state, the oscillators are either locked or 
drifting II 8l 11911 . For given r = lim r ^oor(f), locked oscilla- 
tors have their frequencies in the region |co, | < Jkj?/ (k) and 
thus their phases are locked at sin -1 (C0;(£) / '(/£,■?)) while the 
frequencies of drifting oscillators are in the region |co, | > 
Jkj?/ (k) and thus d§i/dt ^ 0. Drifting oscillators do not con- 
tribute to the order parameter due to the symmetry of the fre- 
quency distribution g(co) = g(— co). Thus the order parameter 
r is computed only in terms of locked oscillators as 



k=l 



2k 



w d&g(&) / dtye">8 



(j) — sin 



-1 



Jkr 



(4) 

The effective field r is obtained by solving the following self- 
consistent equation derived from its definition, 



2tt 



(j) — sin 



«>(*) 
Jkr 



(5) 



(2) To solve Eqs. @ and (|3J, we expand them in powers of r as 



-J>,,-,(^) and^J;^^) , 

(6) 

which are valid for finite N. Here the coefficients are given by 

a% n -\ = (k 2 "- l )u 2 „-i and a 2n -\ = {(k 2n ) / (k))u 2n -i, where 
(k m )=Z lk k m Pd(k), and 



«2n-l 



(-1) 



n— lo — n— 1 



/ 2 («-3/2)! 



«!(«-!)! 



(7) 



from the relation g(co) = (27i)" 1 ' 2 £~ =0 (-co 2 /2)"/n! and 
J% 2 /2 d§ cos 2 $ sin 2 "- 2 (j) = 7t 1 / 2 («-3/2)!/(2n!). Eq. © is 
valid also in the thermodynamic limit unless «2n-i or a% n -\ 
diverge. In case of SF networks with p c i{k) ~ k~^, how- 
ever, (k m ) diverges as N" 1 ^ 1 ^ 1 / (m - k + 1) for m > 
X — 1 1 24] and thus a 2n -\ (« > k/2) and a.2 n -\ (« > 
(k- l)/2) diverge: a 2n -i ~ c^N^/^- 1 /(2n - X) 
for n > k/2 and a 2 „-i ^ c 2 „_iA^ (2n '/ ( ^ 1 ^ 1 /(2n -k+1) 
for n > (k — l)/2 with c 2n -\ and c~2n-\ constants. In 
terms of a scaling variable q = N^'^^Jr/ (k), the terms 
with such diverging coefficients in Eq. © are arranged as 
(J ?/ {k)f- 1 Z„> 1/2 c 2n -iq 2 "- X / (2n -k) = (J ?/ {k)f- 1 A (?) 
for r and (Jr/ (k)) x - 2 ^ n>{x _ l)/2 c 2n ^q 2 "- x + l / (2n - X + 

1) = (Jr/ (k))^~ 2 fx(q) for r. The alternating sign and fast de- 
cay of c 2 n-\ and c 2n -\ from the behavior of u 2n -\ in Eq. Q 
make f\(q) and A(?) finite in the limit q — > °° that corre- 
sponds to iV — > 00 and r finite. Consequently, the following 



equations should be considered for SF networks in the ther- 
modynamic limit, 



n=\ 



-1 



<*>' 



2n-l 



2n-l 



l X-\ 



l X-2 



J 

w 



w 



X-l 



X-2 



, (8) 



where all the coefficients are finite and in particular, a 2 «-i = 

a 2 „-i for n < |A/2J and a' 2n _ l = a 2n -\ for n < [(A- 1)/2J 
with [x\ denoting here the greatest integer that is less than x. 
Notice that d % _ x = f x (°°) and a' x _ 2 = / x (oo) 

The existence of the synchronization transition for X > 3 is 
identified by the emergence of non-zero values of r and 7 only 
when ciyj I (k) > 1 or J > J c with 



C " V K (k 2 ) ■ 



(9) 



Inspecting Eq. (|8jl for 0<A = 7/7 C — 1 <C 1, one finds that 

r ~ ((k) 2 /(k 2 ))r and that Ar ~ |a' 3 |(./ c >7(/fc}) 3 for A, > 5 and 
Ar ~ |fl^_ 2 | (J c F/(k)) 1 - 2 for 3 < A, < 5. Thus it holds for small 
positive A and X > 3 that 



with the order-parameter critical exponent p given by 
P = 



for A, > 5, 
for 3 < X < 5. 



(10) 



(11) 



Notice that in the globally-coupled case where (k 2 ) = (k) 2 = 
(N — l) 2 , the critical point J n is 2 a/2/ 



K ~ 1 .60 and the crit- 
ical exponent P is 1/2 Ill8lll9ll . In Ref. [26], the exponent p 
has been shown numerically to be close to 1 /2 in the synchro- 
nization transition on the Barabasi-Albert network 1 3] that has 
X = 3. 

When 2 < X < 3, (k 2 ) is (aK 3 ^)/^- 1 )) Q and thus the 
critical point/ c is of order N^^^^^ l \ which is zero in the 
thermodynamic limit. It means that r and r are always non- 
zero values for non-zero values of J. Eq. is arranged into 



r~2" 3 / 2 7i 1/2 7r andr: 
leads to 



a x _ 2 \(Jr/(k)) 1 - 2 for small J, which 



1-2 



and 



(12) 



The phase diagram and dynamical states of the oscillators in 
the desynchronized state and the synchronized state are shown 
in Fig. |2 The only contribution to the order parameter in the 
desynchronized state is made by temporary synchronization 
of a few oscillators that could be disconnected. On the con- 
trary, a finite fraction of oscillators have their phases close 
to the average phase 0(f) in the synchronized state. This gi- 
ant cluster of synchronized nodes contains a finite fraction of 
edges as well, implying that the synchrony originates from 
the interaction between them through the SF coupling topol- 
ogy. Furthermore, the continuously-varying critical exponent 



Synchronized 




Desynchronized 

3 4 5 6 7 8 ^ 

FIG. 2: (Color online) Phase diagram and typical dynamical states 
of oscillators on SF networks. The phase boundary in Eq. {5J is 
drawn with Pd{k) = k /t,(X) for k > 1 and ^(A.) the Riemann-zeta 
function. The two network configurations represent dynamical states 
of coupled oscillators that have evolved with J = 1.0 (desynchro- 
nized) and J = 4.0 (synchronized) respectively on SF networks with 
A = 3.6. In both, synchronized nodes, those having their phases (]); 
so close to the average phase 6 that |(|); — 6| < 0.5, are drawn as filled 
boxes inside the circle while the others are as empty boxes at the 
circumference. Edges connecting synchronized nodes are drawn as 
thick lines. 



P according to the degree exponent A demonstrates that the 
variation of the network topology may bring about the change 
of the universality class of the synchronization transition on 
complex networks. The influence of the structural organiza- 
tion on the emergence of dynamical order is more significant 
in case of 2 < A < 3 where J c =0. A deeper understand- 
ing of the nature of the synchronization transition is available 
when we investigate the finite-size scaling behavior of the or- 
der parameter. As we shall see, it illuminates different nature 
of the synchronization transition between SF networks with 
2 < A < 3 and A > 3 as well as allows to confirm our findings, 
Eqs. (BO, (HO), (0, and numerically. 



IV. DERIVATION OF FINITE-SIZE SCALING BEHAVIOR 

In this section, we describe how to compute the order pa- 
rameter in the critical regime and combining it with the results 
in the previous section, present the finite-size scaling behavior 
of the order parameter. 

In finite-size systems, the order parameter r is non-zero 
even with / = so that rj=o ~ N~ l l 2 , as the phases are dis- 
tributed uniformly in [0,27t). With J increasing from / = 
0, the clusters, connected oscillators with nearly the same 
phases, form and evolve to find the largest one among them 
much larger than rj=oN ~ N 1 / 2 , when the order parameter 
r can be approximated by r ~ S/N with S the (ensemble- 
averaged) largest cluster size. The synchronized state is 
marked by the presence of a giant cluster whose size is O (N) 
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FIG. 3: Data collapse of scaled order parameter for (k) = 4, (a) A, = 6.4, (b) X = 3.6, and (c) X = 2.4 according to Eqs. I15i and 1161 . In (a) 
and (b), insets show the crossing of scaled order parameters at J c , which is found to be (a) 1.85(5) and (b) 1.15(5), larger than 1.22 and 0.84, 
respectively, evaluated from Eq. {5J. In the inset of (c), asymptotic behavior of the scaling function is explicitly shown. The solid line has 
slope 1/(3 - A.) = 5/3. 



making r finite. 

To understand the jump to 0(1) of the order parameter in 
the critical regime for X > 3 and in the small J regime for 
2 < X < 3, we should know the largest cluster size S c in those 
regimes. S c is of order N a with 1/2 < a < 1, and we will 
try to extract the exponent a from the asymptotic behavior 
of the cluster size distribution, which can be obtained using 
Eq. (|8), but only for X > 3. For 2 < X < 3, we will see that the 
largest cluster size increases gradually from O (N 1 / 2 ) to (N) 
without any characteristic size between them while J is small. 

Let us define the generating functions 2> (z) = Y,sP( s )z? an d 
<2 (z) — £ s -P(s)z s , where P(s) is the probability that a node be- 
longs to a size-i cluster and P(s) the probability that an edge 
is followed by a size-s cluster. The dependence on the system 
size N is implicitly included. The generating functions are an- 
alytic for finite N and thus one can expand the inverse function 

rP _1 (oj) as z = <P _1 (co) = 1 -L«>i^«(l - co)" around co = 1. 

_ i 

When z = z* N = e 5 with S satisfying S2 <C S <C S and S2 
being the second largest cluster size, the corresponding value 
CO^, = <S (Zff) represents the statistical weight of all clusters but 
the largest one, that is, co^ ~ Ls<s^( s )- Suppose that Eq. l|8} 
is valid in the critical regime. Since f = lim;v^oo(l — COw), 
one can see that the coefficients b n should take such val- 
ues as enable the equation z = 1 — L«>i^n(l — <*>)" with 
z = lirriAr^oo Zff = 1 to be reduced to Eq. HSJi with f replaced 
by 1 — co. Consequently we have z = CO + £" =1 a 2 „_ l [J(l — 
co)/^)] 2 "- 1 +a' x _ 2 [J(l - &)/{k)] x - 2 + ■■■ with co = i (z). 

The last relation informs us that at the critical point where 
fl[7/ (k) = 1, 1 — 2>(z) has the following leading terms: (1 — 
z)'/ 3 for X > 5, (1 -z)V(^-2) for 3 < A. < 5, and J x - 2 (1 - 
z) x ~ 2 for 2 < X < 3, respectively. Eq. (|8j also enables us to 
know that <£ (z) is related to <P (z) via 1 — <P (z) ~ /(l — i (z)) 
and thus the leading terms of 1 — fP (z) are given as (1 — z) 1 ' 3 
for X > 5, (1 - z )V(*.-2) for 3 < A.< 5, and J x ~ l (l -z)^ 2 
for 2 < X < 3, respectively. Here we showed the J de- 
pendence explicitly only for 2 < X < 3. These singulari- 



ties give the asymptotic behaviors of P(s) through the rela- 
tions (l/s\)(d s /dz s )fP(z)\ z=0 =P(s) and {\/s\){d s /dz s ){l - 



T-ll 



s x for large s. Finally the largest cluster size 



S c at the critical point can be obtained using the relation 

L*<s c p ( s ) = 1 - Sc/N jH. The result is 



m for X > 5, 

for 2 < X < 5. 



(13) 



For X > 3, S c in Eq. d 1 3I > is much larger than rj = oN ~ N l l 2 
and we can accept it as true. The oscillators are not ready to 
form the giant cluster until the coupling strength reaches the 
critical value. There emerge clusters of all sizes up to S c at 
the critical point and many of them merge into the giant clus- 
ter with the coupling strength slightly increased. Combining 
Eqs. (II Oi and dl3> . and introducing another exponent /j de- 
fined by S c /N - ArP/^ which is given as 



2 for X > 5, 
'"~ ^ £zi for 3 < A,<5, 



(14) 



we obtain the following finite-size scaling behavior of r 

r = N^ /f,x V(AN l/ ^). (15) 

The scaling function ^(x) is a constant for x <C 1 while be- 
haves as x$ for x — > °°. We remark that o(A^ -1 / 4 ) fluctuation 
of r at J — J c for A, > 5 has also been found in the Kuramoto 
model 1 20] as well as in SW networks l2lll . We plot the scaled 
numerical data for A, = 6.4 and 3.6 in Fig.[3](a) and (b), respec- 
tively, which show good agreement with Eq. (1151 except for 
slight deviations of J c from Eq. 0. Such deviations of the 
critical point have also been observed in the Monte Carlo sim- 
ulation of the Ising model on SF networks 12811 . but are not so 
serious as to give rise to a finite value of J c even for 2 < A. < 3. 
It is known that the critical point can be more precisely pre- 
dicted by studying the system on the Cayley tree with a given 
degree distribution 12911 . 
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When 2 < A. < 3, S c in Eq. Jl 31 is much smaller than rj=oN. 
This means that Eq. is not correct in the critical regime 
because of strong finite-size effects. Instead, the finite-size 
scaling behavior of the order parameter for 2 < A, < 3 can be 
obtained by connecting r ~ J 1 ^ 3- *-) in Eq. (I12> and rj=o ~ 
N~ l l 2 , which leads to 

r = N-?<t>{JN^), (16) 

where 4>(x) is a constant for x <C 1 and <I>(x) ~ ^Vf 3 -^) for 
x — > cxj. This behavior is confirmed by the numerical data with 
A = 2.4inFig.|3](c). 

Distinctive finite-size scaling behaviors shown in Eqs. (1151 
and (I16> represent a significant difference in the emergence 
of dynamical order between SF networks with A > 3 and 
2 < A- < 3. In SF networks with 2 < A < 3, the giant cluster of 
synchrony is formed if only J ^ jsf-^-^)! 2 , A large number 
of hub nodes that are connected with one another do not allow 
separate clusters of various sizes to be developed but force 
most nodes to belong to the unique giant cluster for any non- 
zero coupling strength. It is completely different from the way 
of forming the giant cluster in SF networks with X > 3. Clus- 
ters of various sizes up to (N 1 ^^) in the critical regime 
merge to form the giant cluster as the system escapes from 
the critical regime. The advantage of many real complex sys- 
tems having their degree exponents between 2 and 3 in their 
network structures is thus very obvious. Such highly hetero- 
geneous systems can remain in the dynamically-ordered state 
only with a very small coupling strength, e.g., o(N~^~^/ 2 ) 
in the synchronization phenomena, which is required to main- 
tain the stability of the system. On the other hand, a perturba- 
tion can propagate easily through the hubs as revealed by the 
linear stability analysis of the ordered state 111 5Ul all 711. nec- 



essary to adapt the system to the changes of the environment. 



V. SUMMARY AND DISCUSSION 

We investigated the synchronization transition of coupled 
limit-cycle oscillators on complex networks analytically and 
numerically. We used the mean-field approach to find the 
critical coupling strength and the critical exponent. Further- 
more, we showed how the self-consistent equations derived 
by the mean-field approach can be used also to compute the 
critical fluctuation of the order parameter that completes the 
finite-size scaling analysis. We believe that this method can 
be generally applied to the study of finite-size effects in com- 
plex networks. The synchronization transition turned out to be 
closely related to the formation of the giant cluster of synchro- 
nized oscillators, which crucially depends on the connectivity 
pattern of a given network such that the critical phenomena 
are distinguished according to the degree exponent. Our find- 
ings imply that complex systems take very heterogeneous con- 
nectivity patterns to acquire both stability and flexibility with 
only a small interaction among the elements. While prepar- 
ing this manuscript, we have learned of the work by T. Ichi- 
nomiya [30], which partly overlaps with ours. 
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